## This program generates Figure 9 from "Financial Constraints, Sectoral Heterogeneity,
##  and the Cyclicality of Investment" by Cooper Howes

##### Step 0: Preliminaries #####

## Clear variables and data
rm(list=ls())

## Load packages
library(tidyverse)

##### Step 1: Load data #####

## Load impulse responses from model (created in Dynare file "govspend_model.mod")
irfdata <- t(read.csv("irfs/shock_irfs_govspend.csv",header=F))

## Extract names and relabel
colnames(irfdata) <- irfdata[1,]
rownames(irfdata) <- NULL
irfdata <- irfdata[-1,]
irfdata <- apply(irfdata,2,as.numeric)

## Pick which shocks to show IRFs for
shocknames <- c("e_m","e_a","e_nu","e_sigma","e_xi","e_delta","e_g")

longshocknames <- c("Monetary","Productivity","Labor supply","Risk aversion",
                     "Financial","Depreciation","Gov. spending")

## Plot everything in one figure
windows(width=50,height=50)
par(oma=c(3,3,1,0))
numplots <- length(shocknames)*3
m <- matrix(1:(length(shocknames)*3),ncol=3,byrow=T)
layout(mat = m)

## Generate x-axis labels
qtrs <- 0:16

for(i in shocknames){
  
  irf1 <- irfdata[,paste0("y_",i)] 
  irf2 <- irfdata[,paste0("pi_",i)] 
  irf3 <- irfdata[,paste0("k_",i)] 
  
  y_range1 <- range(irf1)
  y_range2 <- range(irf2)
  y_range3 <- range(irf3)
  
  ## Output panel
  par(mar = c(1.75,3,2,2))
  plot(x=qtrs,y=irf1,type="l",col="blue",lwd=3,ylim=y_range1,ann=FALSE,cex.main=1.25,
       axes=F,frame=T)
  axis(1,cex.axis=1.5)
  axis(2,at=pretty(y_range1,n=5),lab=paste0(pretty(y_range1,n=5)*100, "%"),cex.axis=1.25,las=TRUE)
  abline(a=0,b=0,col="black",lwd=1)
  mtext(text=longshocknames[which(shocknames==i)],side=2,col="blue",line=4.5,cex=0.8,font=2)
  if(i=="e_m"){title(main="Total output", col.main="blue", font.main=2,cex.main=1.25)}
  
  ## Inflation panel
  par(mar = c(1.75,3,2,2))
  plot(x=qtrs,y=irf2,type="l",col="blue",lwd=3,ylim=y_range2,ann=FALSE,
       axes=F,frame=T)
  axis(1,cex.axis=1.5)
  axis(2,at=pretty(y_range2,n=5),lab=paste0(pretty(y_range2,n=5)*100, "%"),las=TRUE,cex.axis=1.25)
  abline(a=0,b=0,col="black",lwd=1)
  if(i=="e_m"){title(main="Inflation", col.main="blue", font.main=2,cex.main=1.25)}
  if(i=="e_g"){mtext(text="Quarters after shock",side=1,line=3)}
  
  ## Capital stock panel
  par(mar = c(1.75,3,2,2))
  plot(x=qtrs,y=irf3,type="l",col="blue",lwd=3,ylim=y_range3,ann=FALSE,
       axes=F,frame=T)
  axis(1,cex.axis=1.5)
  axis(2,at=pretty(y_range3,n=5),lab=paste0(pretty(y_range3,n=5)*100, "%"),las=TRUE,cex.axis=1.25)
  abline(a=0,b=0,col="black",lwd=1)
  if(i=="e_m"){title(main="Capital stock", col.main="blue", font.main=2,cex.main=1.25)}
  
  
}

